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An ill-posed inverse problem of autoconvolution type is investigated. This inverse problem oc- 
curs in nonlinear optics in the context of ultrashort laser pulse characterization. The novelty 
of the mathematical model consists in a physically required extension of the deautoconvolu- 
tion problem beyond the classical case usually discussed in literature: (i) For measurements 
of ultrashort laser pulses with the self-diffraction SPIDER method, a stable approximate 
solution of an autocovolution equation with a complex-valued kernel function is needed, 
(ii) The considered scenario requires complex functions both, in the solution as well as in 
the right-hand side of the integral equation. Since, however, noisy data are available not only 
for amplitude and phase functions of the right-hand side, but also for the amplitude of the 
solution, the stable approximate reconstruction of the associated smooth phase function rep- 
resents the main goal of the paper. An iterative regularization approach will be described that 
is specifically adapted to the physical situation in pulse characterization, using a non-standard 
stopping rule for the iteration process of computing regularized solutions. The opportunities 
and limitations of regularized solutions obtained by our approach are illustrated by means 
of several case studies for synthetic noisy data and physically realistic complex- valued kernel 
functions. Based on an example with focus on amplitude perturbations, we show that the 
autoconvolution equation is locally ill-posed everywhere. To date, the analytical treatment 
of the impact of noisy data on phase perturbations remains an open question. However, we 
show its influence with the help of numerical experiments. Moreover, we formulate asser- 
tions on the non-uniqueness of the complex-valued autoconvolution problem, at least for the 
simplified case of a constant kernel. The presented results and figures associated with case 
studies illustrate the ill-poscdness phenomena also for the case of non-trivial complex kernel 
functions. 



1. Introduction 



In the early 1990s motivated by applications from spectroscopy (cf. an d stochas- 
tics (cf., e.g., [2, p. 74]) contributions to the deeper mathematical and numerical 
analysis of deautoconvolution problems as a class of inverse problems in spaces of 
continuous or quadratically integrable real functions were made. Precisely, deau- 
toconvolution problems under consideration were mostly aimed at finding non- 
negative functions x with compact support supp(x) C R from its autoconvolution 
x * x. After transformation to the unit interval supp(x) C [0, 1] they consist in 
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finding solutions of the autoconvolution equation 



x(s - q)x(q)dq = y(s) , (1) 



where the support of y belongs to the interval [0,2]. Since the autoconvolution 
operator x i— > x * x is nonlinear and 'smoothing', the deautoconvolution problem 
is ill-posed in the sense that for given y the solutions x need not be uniquely 
determined and mainly small perturbations in the right-hand side y caused by noisy 
data may lead to arbitrarily large errors in the solution. To overcome the negative 
consequences of ill-posedness up to some extent some kind of regularization is 
required. Regularization techniques allow us to find stable approximate solutions of 
equation ([JJ based on auxiliary problems. For data y on the subinterval < s < 1 
the paper |3j has analyzed the situation of equation ([I]) and the application of 
Tikhonov's regularization method including its convergence properties, whereas in 
(H the situation of data y on the whole interval < s < 2 has been studied. 
Alternative regularization methods applied to equation ([1]) and specific numerical 
approaches were also discussed in @-ll]. 



Recently, the research group 'Solid State Light Sources' of the Max Born In- 
stitute for Nonlinear Optics and Short Pulse Spectroscopy, Berlin, hit on the au- 
toconvolution problem in the context of a new approach in ultrashort laser pulse 
characterization called Self-Diffraction SPIDER or short, SD-SPIDER (cf. [13]). 
For phase reconstruction as an auxiliary problem the solution of equation (JTJ) is 
needed, but for complex functions x : [0, 1] C R — > C to be determined from com- 
plex observations y : [0, 2] C IR — > C. To our knowledge, a thorough analysis of the 
complex case in deautoconvolution is still missing in the literature, in particular 
as the ill-posedness phenomenon arises in the complex case, too. Moreover, the 
occurrence of a device-related kernel function k : [0, 2] x [0, 1] — > C involved in 
the mapping x i— > y, which is non-trivial in the sense that k ^ 1, constitutes a 
challenging additional difficulty in connection with this inverse problem behind SD 
SPIDER. So as a part of the SD SPIDER approach it would be necessary to solve 
(after transformation of the variables to the unit interval) the equation 

s 

J k(s,q)x(s-q)x(q)dq = y(s), < s < 2, (2) 


in a stable approximate manner when only noisy data of y are given. To simplify 
the notation we write integrals like on the left-hand side of ([2]) even if a function 
in the integrand is not defined there, as this is the case for x(s — q) if s — q > 1. 
Then we set the corresponding function values as zero and avoid to distinguish the 
integral representations for < s < 1 and 1 < s < 2. 

The equation ([2]) is a complex-valued and kernel-based generalization of equation 
([T]) with solution 



and right-hand side 



x{q) = A(q) exp[i<p(q)], < q < 1, (3) 



y(s) = B(s) exp[#(s)], < s < 2. (4) 



We consider in this paper two different aspects of solving the integral equation ([2]) 
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as a mathematical model for inverse problems. The general inverse problem (a) 
consists in finding the complex function x in ([2]) from noisy data of the complex 
function y and for a given complex- valued kernel k. Alternatively, the SD-SPIDER- 
motivated specification (b) lies in finding the phase function ip in ([3]) from noisy 
data of y and given k when additional observations A of the amplitude function 
A are available. Below we will emphasize that the identifiability of the phase <p 
requires measurement data of both the phase function tp as well as its correspond- 
ing amplitude function B in In order to find the phase function tp in problem 
(b) under some smoothness assumptions we will suggest a Tikhonov regulariza- 



tion approach (cf., e.g., [13|, LLJ]) combined with a specifically adapted strategy for 



choosing the regularization parameter a > 0. In this context, iterative procedures 
of Levenberg-Marquardt type (cf., e.g., [IH]) are required to compute acceptable ap- 
proximations of the Tikhonov-regularized solutions in an efficient manner. Exploit- 
ing the additional knowledge of noisy data for the amplitude function A allowed 
us to construct an adapted stopping rule for the developed iterative regularization 
approach. Moreover, we illustrate the quite acceptable work of this approach by a 
numerical case study. 

Currently, not all questions concerning (a) and (b) can be answered by the au- 
thors. Therefore, it makes sense to show the local ill-posedness of problem (a) by 
an analytic example and to illustrate ill-posedness phenomena of (b) by numerical 
case studies. At the moment the analytical treatment of the impact of noisy data 
on phase perturbations is an open question. Moreover, it seems to be reasonable to 
present a uniqueness assertion on the inverse problem (a) in the case k = 1 since 
assertions for general complex kernels are missing to our knowledge. On the other 
hand, all numerical case studies for problem (b) were performed with physically 
relevant kernel functions k. 

At this point we should note that the problem of finding a function ([3]) solving the 
equation ([2]) is a generalization of the problem of recovery of a compactly supported 
and complex-valued function from the modulus of its Fourier transform. This so- 
called phase retrieval problem and its applications, for example in optics, electron 
microscopy and astronomy, were intensively studied in the literature based on the 
seminal paper [lfj|, and we refer to the review paper [H], to [H, Section 1.2], and 
references therein. In particular, rigorous uniqueness results were proven for the 
one- and two-dimensional phase retrieval problem (cf., e.g., 16, 3-2lt]). Moreover, 



numerical methods were proposed in [171 . 



2fl]. 



The paper is organized as follows: in Section 2 we briefly review the role of 
deconvolution and decorrelation in the characterization of laser pulses. Then we 
analyze in some more detail the previously unsolved deconvolution problem in self- 
diffraction SPIDER, i.e., a specific variant of laser pulse characterization methods. 
Leading from the notation commonly used in the physical literature to the one em- 
ployed in mathematics, we reformulate the abstract mathematical problem behind 
and investigate its ill-posedness in Section 3. Subsequently, in Section 4, we derive 
consequences of the Titchmarsh convolution theorem for a constant kernel function 
k and present an adapted regularization approach in Section 5. Finally, we consider 
several case studies based on synthetic noisy data and illustrate the convergence 
behavior of the iteration procedure, which requires a non-standard stopping rule. 
After a summary of the autoconvolution problem and its specific regularization 
approach, we conclude with a brief outlook at the relevance of the findings in the 
physical sciences. 
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Figure 1.: (a, left) Illustration of the pulse to be characterized. The shortest gener- 
ated laser pulses comprise only few oscillation of the electrical field £(t). The aim 
of pulse characterization is to retrieve the temporal evolution of the envelope of the 
electric field A(t) or of the intensity I(t) = |^4(t)| 2 . (b, right) The Fourier domain 
representation of the pulse is often utilized within the process of pulse retrieval. 
Therefore, it is decomposed into its amplitude A{uj) and the spectral phase ip(uj). 



2. Physical background 

2.1. The evolution of ultrashort pulse characterization 

Ultrashort laser pulses constitute the shortest man-made controllable events, with 
demonstrated pulse durations of 4 femtoseconds (4 x 10" 15 seconds) and below 
224241] . Using wave leng th conversion techniques, pulses as short as 8 x lCP 17 s 
have been produced [25l. l26|. which is in the range of the fastest transient events in 
atoms and molecules. Quite remarkably, these generated pulses approach a limiting 
width of a single optical cycle of the underlying electric field carrier, as illustrated 
in Figure [H Nevertheless, these remarkable technological achievements also cause 
a serious dilemma for their accurate characterization and measurement. Temporal 
resolution of a dynamical process is always limited by the gate time of the sampling 
process, similar as in photography, where temporal resolution is dictated by fastest 
available shutter speed. As there simply exists no gate process faster than the 
duration of the shortest laser pulses, the history of laser pulse characterization 
and ultrafast spectroscopy has also always been a history of deconvolution and 
decor relation. 

From a physicist's point of view, a practical and easy-to-implement remedy for 
the fundamental dilemma of pulse characterization is the use of identical replicas 
of the input pulse as the pulse under investigation and for the gate function (2?J • 
For implementation of the temporal gate, one has to form the product of the two 
functions, which is conveniently done by employing a nonlinear optical process. In 
the simplest case of sum frequency generation, this process generates the product 
of two intensity envelopes X\{t) and ^(t), with one of these envelopes serving as 
the gate event. Varying the delay between these temporally dependent signals then 
allows recording their cross-correlation 



C(r) oc / lx(t)l 2 (t-r)dt, 



(5) 



where we employed the usual sign convention for the delay r in optics and the 
integration over time t resembles the integration that is carried out by the slow 
detector recording the correlation signal. Reconstruction of Ii(t) with known corre- 
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lation function C(t) and gate function 2a(t) constitutes an inverse problem, which 



is frequently solvable with satisfactory precision, see, e.g., [28|] for a practical ex- 
ample. It is important to understand that experimental noise mandates the gate 
function to be kept as simply structured as possible. Ideally, it should be chosen a 
single-maximum function with smallest possible temporal duration. 

Using identical functions Z{t) = T\(t) = Z2(t), Eq. ([5]) describes an autocorre- 
lation, and the decorrelation problem turns out to be only ambiguously solvable 
[1^,[30]. This can be easily seen by use of the Fourier convolution theorem. In the 
spectral domain, the correlation is written as a product 

C(u) oc I(u)I*(u) = \I{oj)\ 2 , (6) 

where u) is the optical angular frequency and I(oj) is the spectral intensity. The 
Fourier transform f(uj) of a function J-(t) is given by 

00 

f(u)= J F(t) exp (-tut) dt. (7) 



In Eq. ([6]), the phase of the complex- valued I(oS) does not affect the correlation sig- 
nal C(oj). This means that all information on temporal pulse asymmetry in the time 
domain does not enter into the autocorrelation signal. Simplistic reconstruction, 
trying to revert the absolute square in ([6]) by the principal square root, therefore 
always yields a symmetric reconstructed Z(t) . Consequently, decorrelation does not 
unambiguously work unless additional information on pulse asymmetry is obtained. 



Several ways have been investigated to resolve this issue. Early attempts [3l|, [32|] 



relied on the availability of additional spectral information I(uj), further attempts 



included more complex types of autocorrelations 33J]. It could be shown, however, 
that even the tiniest amounts of measurement noise render all these concepts un- 
practical [33 ]. yet again thwarting unambiguous reconstruction of I(t). In all these 
types of decorrelation problems, therefore, thorough analysis typically pointed out 
severe limitations rather than offering practical solutions. 

This unsatisfactory situation changed dramatically with methods like frequency- 
resolved optical gating (FROG [H]) and the sonogram technique [3(J. In simple 
words, the idea of these method is to replace the separate acquisition of a spec- 
trum and an autocorrelation by measurement of a two-dimensional spectrogram- 
like function S(t, uj) with simultaneous functional dependence on delay and fre- 
quency. In the simplest case of second-harmonic FROG, 

00 

S(t,u)(x J l(t)l(t - t) exp (-iwi) dt. (8) 



In the experiment, the acquisition of FROG traces only requires replacing the 
spectrally integrating detector by a spectrally resolving one and can be used with all 



existing autocorrelation geometries. It can be shown [35[ that FROG resolves most 
but not all issues of decorrelation. There is a one-to-one correspondence between 
a given FROG trace S(t, u) and the intensity envelope I(t), provided that the 
pulse does not consist of temporally or spectrally separated segments with intensity 
levels approaching the experimental noise floor in between. A detailed discussion 



of remaining ambiguities can be found in [33, |38( . One additional ambiguity of a 



second-order autocorrelation- type FROG method is time reversal, i.e., the method 
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cannot distinguish between Z(t) and I(—t). Nevertheless, there are alternative 
methods relying on third-order nonlinearities [35| and other fixes [39I ] to overcome 
this issue. Retrieving I(t) from a known S(t, co) is yet again an inverse problem, 
and considerable effort went into reliable and numerically efficient algorithms for 
retrieving the pulse shape from a FROG trace (3BI . 4ol- 44 ] . 



As an alternative to FROG, spectral phase interferometry for direct electric- 
field reconstruction (SPIDER) emerged 45, [4f|. This method is conceptually dif- 
ferent from all previous approaches and normally does not require any kind of 
deconvolution. Being based on spectral interferometry [47I |4^ |. SPIDER directly 
measures the spectral phase (f(co) of the electric field E(lo) = A(u)exp[iip(uj)], 
where both the spectrum A{oj) and the phase (p(ui) are real-valued functions. In 
the time domain, the field can again be separated into amplitude and phase via 
£(t) = A(t) exp [i<fi(t)]. The electric field relates to the intensity via Z(t) = |^4(t)| 2 . 
SPIDER requires the generation of two replicas of the input pulse X(i), which are 
processed with an ancillary pulse X a (i) in a second-order nonlinear crystal. This 
ancillary pulse can also be derived from the same input field, i.e., the method is 
self-referenced. To this end, a fraction of the input pulse is sent through a dis- 
persive medium, e.g., a glass block. The dispersion of the glass block gives rise 
to a non-vanishing chirp = dcf)(t)/dt, i.e., the pulse's carrier frequency varies 
with time. Mixing the field of the ancillary pulse <? a with the two replicas con- 
tained in £ r causes a different frequency shift for each replica, which is referred to 
as spectral shear. The spectrum of the nonlinear mixing product, i.e., the Fourier 
transform of £ r (t)£ a (t) can then be analyzed with the methods described in 4^, 50] 
to reconstruct the spectral phase <p{u) of the original input pulse. With additional 
knowledge of the amplitude A(oS), a simple Fourier transform then suffices for 
complete reconstruction of I(t) in the time-domain. 

While traditional SPIDER methods all rely on second-order nonlinearities, it 
seems appealing to extend this proven method to third-order nonlinearities as they 
are frequently used for autocorrelation and FROG. These higher-order methods 
allow generating a FROG or SPIDER signal that is spectrally collocated with 
the generating wave, which alleviates pulse characterization in the ultraviolet 
Moreover, third-order SPIDER has recently been demonstrated using a monolithic 
waveguide device [5l[ , bringing us closer to the dream of an all-optical oscilloscope 
on a single optoelectronic chip. As third-order optical nonlinearities appear in any 
kind of dielectric material and not only in non-isotropic crystalline media, their use 
is essential for optically integrated characterization methods. Unfortunately FROG, 
autocorrelation and also a SD-based spectral-interference pulse characterization 
technique 54| involve mechanical scanning of optical delays, which rules them out 
for integrated optical devices. At the current state of the art, this only leaves 
SPIDER for such applications. 

A third-order SPIDER, however, turns out to be difficult to implement 
Given that there is now the product formation of three waves involved, which in- 
teract in the nonlinear mixing process, one can either generate two waves from the 
ancillary pulse and one from the two replicas, or, vice versa, two replica waves £ T (t) 
and only one ancilla £ a . As the chirped and temporally stretched ancilla typically 
displays a much lower peak intensity than the two replicas, the latter approach 
with two replica waves and one ancilla therefore constitutes the more favorable 
physical scenario, resulting in higher conversion efficiencies. Using, for example, a 
SD geometry the nonlinear mixing process forms a signal £sd(*) °t £ r 2 (^)£a(*)- 
This duplicates the scenario of second-order SPIDER with two notable differences. 
The replicas enter squared £^{t), and the ancilla enters complex conjugated. While 
the latter essentially only corresponds to a sign reversal of the squared repli- 
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cas cause SPIDER to measure the spectral phase ( / 9 conv (w) of the convolved signal 
£conv{t) oc £ r 2 (i) rather than the spectral phase of the pulse itself. Hence, a decon- 
volution is required in order to gain insight about the shape of the electric field 
£(t) or of the intensity envelope T{t). 



2.2. The deconvolution problem in SD SPIDER 

The derivation of the deconvolution task formulated in Eq. (2) - (4) requires a more 
detailed analysis of the generation process of the SD wave. In the last paragraph, 
a simplifying time-domain-based formulation of the nonlinear signals generated for 
pulse characterization has been employed in order to work out the basic differences 
between the different characterization methods. A more general formulation is ob- 
tained in the Fourier domain and by the use of the nonlinear wave equation for the 



involved electro- magnetic waves [53]] . From this approach one can derive under the 
slowly varying envelope approximation that the strength of the SD signal is given 
by 



E SD (L,oj) = i^ j P n ,(oj)e~^^ z dz, (9) 

where P n i(z,ui) is the nonlinear polarization of the bound electrons inside the 
materials employed for generation of the SD signal. The vacuum permeability /iq, 
the speed of light c and the linear refractive index n = n{u) are physical constants. 
z denotes propagation direction of the electro-magnetic waves inside the material 
that is used to generate the SPIDER signal and L denotes the material's length. 
u) is the frequency of the self-diffraction signal and ksd(^) is the wave number of 
this signal. The constants in Eq. ([9]) and the integral over the interaction length 
will contribute to the kernel k of the autoconvolution problem Eq. ([2]). For further 
analysis the nonlinear polarization P n \(uj) has to be considered, which is in general 
given by 

roo roc 

-PniM = / / x^\u,-u) a ,u)i,u)2)ElE I (co 1 )E t (u) 2 )6(u) + u) a ,-u)i - lo 2 ) dw 2 duii . 
Jo Jo 

(10) 

Here, uji, u 2 are the frequency components of the two replica waves E r , and a> a 
denotes the single frequency component of the ancilla E a . The values of x 
and 5{u + w a — oj\ — 0J2) express physical constraints and relations of the light- 
matter-interaction: x^ 3 \ the third-order nonlinear susceptibility, denotes the in- 
teraction strength between the field of light and the material in use. The factor 
S(u) + w a — uj\ —L02) ensures energy conservation in the expression above. Inserting 
of Eq. (|10p in Eq. ([9]) and simplification of the result yields 

-EsdM = / x (3) (w,-w a ,wi,w + w a -a;i)M(a;,-u; a ,a;i,a;2) 

J ° (ID 
E*E r (u;i)E r (uj + w a — uji) dwi . 



Here, M collects the constants and the result of the z-integration of the exponential 
term in Eq. describing the phase matching conditions between all waves involved 
in this four-wave-mixing process: The frequency dependent values of an d M 
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are known from experiments and merge to the kernel 

iTKwi.Wa) = ^x^H^-Wa^^+Wa-wiJ^e^-^'TrsincCA^L^) (12) 

where Ak = (Ak^, Ak^, Ak z )(w, -w a ,wi,w+w a -wi) is the phase-mismatch vector 
and (£, 77, z) denotes spatial coordinates. Now ([lip can be written as 

-EsdM = / if(w,wi,w a )£' r (wi)£; r (a; + a; a -wi)(ia)i. (13) 
J o 

We note that the wave i£ r has a spectrum of finite width, i.e., there exist frequencies 
< uf m < wf 1 ^ < oo such that £ r (wi) = for wi ^ [^f^wf 8 *]. Consequently 
.Bsd is a l so compactly supported and E'sd(w) = for uj ^ [2a;f in — w a , 2o;!j nax — u; a ]. 
Substituting 

w + w a = suf 3 * for s G [0,2], 
Wl = go;f iax for g G [0,1] 



in Eq. (fT3]l and defining 

a;(g) 



^ r (^r x ), (14) 
^ SD ( SW r x -w a ), 



allows us to trasnform the limits of the integral such that we arrive at the ab- 
stract mathematical model equation which will be studied in the following. 
By designing the SPIDER apparatus in such a way that one of the replica beams 
can be blocked, the SD signal's amplitude B(s) and its phase ip(s) in ([!]) can be 
measured with the same device such that (noisy) data of the complex function 
y( s )i < s < 2, are available. An independent measurement of the spectral in- 
tensity I(oj) = \A(lo)\ 2 in front of the SPIDER apparatus allows us to compute 
approximations A(q) of A(q) corresponding to the incident pulse. Since for the 
retrieval of the pulse shape, the complex function x(q), < q < 1 has to be 
Fourier transformed, the completely unknown phase function <p(q) or its derivative 
GD(q) = <p'(q) called group delay remains to be determined. For simplicity, in 
the following we neglect the uncertainty of kernel data and suppose to know in a 
precise manner the continuous complex function k(s,q), (s,q) G [0,2] x [0,1]. We 
denote by 

&max := max \k{s,q)\ 

(a,g)e[0,2]x[0,l] 

the maximum of its modulus. 



3. The abstract mathematical model and its ill-posedness 

The physical inverse problem under consideration described in Section [2] requires 
the solution of the generalized autoconvolution equation ([2]). As outlined in Sec- 
tion [H the two aspects (a) and (b) are under consideration. As is well-known the 
comprehension of additional information about expected solutions and the retrieval 
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of data plays an important role for the stable approximate solution of inverse prob- 
lems. Therefore, physicists are preferably interested in aspect (b) aimed at finding 
the continuously differentiable phase function cp : [0, 1] — > R in ([3]), given the kernel 
k and noisy observational data of y and A. We note that the additional knowledge 
of an estimate A(q), < q < 1, of the modulus function \x(q)\ = A(q), < q < 1, 
in ([3]) will play a prominent role for choosing the regularization parameter in the 
process of constructing stable approximate solutions to (|15p in Section [5j Neverthe- 
less, it will give some insight into the problem structure to focus in this and in the 
subsequent section on aspect (a), where ([2]) is considered as a nonlinear operator 
equation 



F(x) = y, (15) 



formulated in appropriate abstract function spaces. 
Taking into account that the forward operator attains the form 



[F(x)}(s):= k(s,q)x(s-q)x(q)dq, < s < 2, (16) 



we are searching for the function x : [0, 1] — > C from noisy data y s G L 2 (0, 2) of y 
that satisfy the deterministic noise model 

\\y s -y|U=(o,2) < S, (17) 

with noise level <5 > 0. Our focus is on the Hilbert space situation F : L 2 (0, 1) — > 
£2(0,2), where F is mapping between Hilbert spaces of square-integrable complex 
functions. For x : [0, 1] C, y G [0, 2] C, and k : [0, 2] x [0, 1] C formula (US} 
is an abbreviation of the form 

s 

J k(s, q) x(s — q) x(q) dq, < s < 1, 
[F(x)](s) := { \ 

J k(s, q) x(s — q) x(q) dq, 1 < s < 2, 

s-l 



which we always use for simplicity assuming that x, y can be extended to IR and 
k to R 2 as zero outside of the original domains. In this context, we remember the 

1 /2 

structure of the norm ||z||l2( 0)£1 ) = (/ Q a \z{q)\ 2 dq) for elements z G L 2 , (0,a). If we 
consider square-integrable real functions z G L 2 (0,a), then the norm is the same, 
but \z(q)\ denotes the modulus of the real value z(q). 

The ill-posedness phenomena of non-uniqueness and instability well-established 
for the equation CD (cf. @, i]) also occur when the complex-valued function x : 
[0, 1] — > C is determined from Eq. ([2]) and also in case that only ip has to be found. 
Together with x also the function — x having the same modulus \x\ satisfies ([2]). 
Hence, if a continuously differentiable function ip : [0, 1] — > R as a part of x in ([3]) 
can be chosen such that equation ([2]) is satisfied for given functions A and y, then 
ip + 7r also solves the equation. On the other hand, together with tp also ip + 2tt and 
hence all tp + run, m G Z, solve the equation. This type of non-uniqueness cannot 
be neglected, but is of inferior significance, since the corresponding group delays 
(p' are uniform for all integers m. 

Instability is a more important difficulty occurring in all linear and nonlinear in- 
finite dimensional least-squares problems which are aimed at solving inverse prob- 
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lems with smoothing forward operators. Therefore, extreme care must be exer- 
cised when discretizing an infinite dimensional least-squares problem since the 
finite-dimensional approximating least-squares solutions may not converge (cf., 
e.g., [55j j) or, even worse, they may diverge from the true solution with arbi- 
trarily high speed (cf., e.g., [56(]). For the nonlinear operator F from (| 16|) we 
unfortunately have that for every xq £ F\ (0, 1) there exist, in case of arbitrar- 
ily small radii r > 0, sequences {x n }^ =1 C L 2 (0, 1) n B t {xq) in a neighborhood 
B r (x ) := {x £ L\ (0, 1) : \x - £o|Ug(o,l) < r } of ^0 with 

\\F(x n ) - F(x )\\li(o,2) 0, but \\x n - scoIIl^o.i) A as ^-^oo. (18) 

We call this local ill-posedness at the point xq and it has the consequence that 
a solution of equation ([2]) cannot be ap pro ximated arbitrarily good even if the 
noise level of the data tends to zero (cf. [57H59I]). For injective operators F local 
ill-posedness indicates that F is not continuously invertible, and this mostly results 
from compactness of the forward operator F. However, in 0] it was shown that 
the autoconvolution operator from equation ([I]) is locally ill-posed everywhere in 
L 2 (0, 1), but fails to be compact. So compactness also cannot be expected as an 
intrinsic property of F in the complex case, but following an idea from 0], which 
was extended in [gO], we can nevertheless prove local ill-posedness ([18]) of F from 
([2]) everywhere in L 2 (0, 1) by Example 13.11 

Example 3.1 For any radius r > the positive function ^^(q) '■= r ^ 1 q ^ 2 ^ , 
< q < 1, possesses for all < j3 < | the properties 

^a 2 (0,l)c^(0,l) and 11^11^(0,1) =r- 
Then the equalities 

[*p**p](s) = r 2 (l-2/3)s 1 " 2/3 f (l-u)- p u- p du = r 2 (l - 2/3)s 1 " 2/3 J B(l - 0, 1 - 0) 



with Euler's beta function £?(•,•), satisfying for < f3 < 1/2 the condition 
< B(l - (3, 1 - /3) < 7r, yield the estimate 



11*0 * *p\\LUo,2) < max [*/9 * */?K s ) < >/2r 2 (l " 2/3)vr2 1 -^ -> as /3 -> - 
By setting x n := xq + we have for every xq G L 2 (0, 1) and all < s < 2 



|[F(x n )](s)- [F(x )](s)\ < I (k(s,q) + k(s,s-q))x (s-q)*L_i(q)dq 

s 

+ J k(s,q)^i_i(s - q)^i_i(q)dq 


s 

< k max {2 / \x (s - (q)dq + * *i_i](s)} ->■ as n oo, 

o 

because is weakly convergent in L 2 (0, 1) to the zero function as n — >■ oo. 
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Since moreover — [F(a;o)](s)| is bounded from above by a constant that 

does not depend on s, Lebesgue's dominated convergence theorem implies that 
\\F(x n )-F(xo)\\ L 2( 0i 2) -> as n -> oo and hence, with ||x n -x ||L=(o,i) = r > 0, ((TSJ 
is valid. Consequently, we have local ill-posedness for all functions xq G L^,(0, 1). 

Note that with respect to x(q) = A(q) exp[i(p(q)] the local ill-posedness men- 
tioned in Example 13.11 refers to exploding amplitudes A(q) for small q as a con- 
sequence of the fact that ^p(q) has a weak pole at q = 0. If A(q) is given as a 
continuous function, then the ill-posedness with respect to the remaining unknown 
continuous phase function <p is less obvious, but case studies show that instability 
also occurs in the sense that clearly distinguished phase functions (p can lead to 
nearly the same complex-valued function [F(x)](s) = B(s) exp[iip(s)]. Such stud- 
ies also prove that <p is not identifiable alone from the phase function ip. Both 
functions, amplitude B and phase ip, are required to recover <p when A is known. 



4. Titchmarsh convolution theorem and its consequence 

For any function x G L^(0, 1) the elements F(x) and F(—x) according to the 
operator F from (|16p coincide and belong to L^(0, 2). However, it is of interest 
whether this is the only ambiguity of solutions to equation ([2]). For the special 
case k = 1 a positive answer can be given by the Titchmarsh convolution theorem 
(cf. [oil ]) which we formulate as a lemma: 

Lemma 4.1: Let f,g£ L^.(R) with supp(/) C [0, oo), supp(g) C [0, oo), and let 
for some constant a > 

s 

[f * g](s) := J f(s — q)g(q)dq = for almost all s€[0,a]. 
o 

Then there are non-negative constants a\ and a<i such that a\ + ai > a and 
f(q) = for almost all t £ [0, oi], g(q) = for almost all q G [0,02]. 

This lemma allows us to prove the following theorem: 

Theorem 4.2: If for given y G L^(0, 2) the function x G L^(0, 1) solves the 
equation (0) with k = 1, then x and —x are the only solutions of this equation. 

Proof: Let x G ^(0,1) and x + A G L 2 c (0,l) solve for all < s < 2 the 
equation (pQ). Then we have [(x + A) * (x + A) - x * x](s) = [A * (2x + A)](s) = 
for almost all s G [0, 2]. If we set 

<7A := sup{g > : A(r) = for almost all r G [0, q]}, 

then x and x + A are different elements of I? c (0, 1) if and only if a\ := q/\ < 1. For 
that case, Lemma 14.11 ensures with 02 > 1 that [2x + A](g) = and [x + A](g) = 
[— x](q r ) for almost all q G [0, 1]. Thus —x is the second solution besides x and other 
solutions can be excluded by this proof. □ 

It seems to be an open problem under what conditions imposed on k ^ 1 the 
result of Theorem 14.21 of having just a twofold solution of ([2]) can be formulated 
and proven. 
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5. An adapted regularization approach 

Now we return to the SD-SPIDER-motivated aspect (b), where additional data 
A are available. To find stable approximate solutions for the operator equation 
(I15|) with the nonlinear operator F : L^(0, 1) — > L^(0,2) from (]16p . we can ex- 
ploit the nonlinear Tikhonov regularization (see, e.g., [13|, Chapter 10]), where the 
regularized solutions x^ a £ L 2 (0,1) are minimizers of 

\\F(x) - ^11^2(0 2) + aSl(x) -> min, subject to x £ D(Q) C L 2 C (0, 1), (19) 

with a regularization parameter a > and a stabilizing functional 0, : D(£l) — > 
[0,oo) (cf. (l8l . Chapter 4]) having the domain D(£l). To prefer smooth solutions x 
the penalty functional is frequently set as 

ft(x) := \\x — x\\ 2 L 2( ^ and i}(x) := \\Lx\\ 2 L 2^ ^, respectively, 

where x G ^(0, 1) is a reference element and ||x — ^H^g ^ attains small values 

if x is close to x, and on the other hand L : D(L) C L^(0, 1) — >• -^(0, 1) denotes 
a densely defined differential operator such that ||£x|| 2 attains small values if for 
example first or second derivatives of x are 'small'. For the calculations we chose 
the second derivative Lx = -§px. The specific measurement situation of our inverse 
problem in ultrashort laser pulse characterization, where only the phase ip in x is to 
be determined, whereas an observation A = \x\ of the amplitude can be observed, 
allows us define the problem specific rule for choosing the regularization parameter, 
which in continuous formulation reads as 



dq for all a > 0, 
(20) 

and to use as the adapted approximate solution to Eq. (fT5|) . Owing to that 
additional data information A the common search of a > based on heuristic 
rules like the quasi-optimality rule, for example successfully applied in [(32|], can be 
completely avoided. Although analytic properties of (j20j) . especially existence of a 
minimizer a*, are hard to verify for a search over all a > 0, the situation simplifies 
in the practically relevant case of a set of discrete values cti, ■ ■ ■ , ckjv} with 
n 6 N fixed. A minimizer of (120p over only the ay, j = 1,...,N, exists. In the 
unlikely case that a* is not unique, we may pick one of the minimizers. Since the 
total amount of computational work for obtai ning xl from (1191) and (120p is rather 
high, iterative regularization procedures (cf. [la, |63( for an overview) can yield 
alternatives to the Tikhonov regularization with reduced computational expenses. 
Our focus is on a variant of the Levenberg-Marquardt method, which is a Newton- 
type method for nonlinear least-s qua res problems. For the mathematical theory of 
this method see the recent paper 164- Here, we consider the iteration process 

8 . I ( rpf f '§ \* rpf f sy.8 \ | _ r * r \ 771/ ( 8 \# 77 / ™5 



a* = a*(y d ,A): / dg < / \ x d a {q)\ - A(q) 



with appropriate relaxation factors 7 > and a regularization parameter a > 0, 
aimed at minimizing the linearized functional 

11/ - F(X(l)) ~ F '( X (l))( X ~ X (l))\\ll(0,2) + a \\ L ( x - X (l))\\ll(0,2) 



+ 7 [F'{xLyF'(xL) + aL*L) F'(x 5 (l) f(y 5 - F{x*% (21) 
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and terminated for I = Z* early enough according to some stopping rule. The last 
iterate x s a := x^ >. acts as regularized solution. In contrast to the classical version of 
the Levenberg-Marquardt method we do not diminish the regularization parameter 
a with growing I £ U, but keep it constant. Under all such regularized solutions 
we select by the parameter choice rule ([20]) and use it as approximate solution 
to Eq. (|15p . Numerical experiments proved for our problem that the squared de- 
viations of absolute values occurring in (|20|) are also helpful for any fixed a > 
to define the index Z* for stopping the iteration process. Precisely, we always ob- 

1 2 

served that the values of discretized versions of the integral J — A(q) dq 

o 

decrease with growing Z = 1, 2, ... in the initial part of the iteration up to some turn 
around point with iteration number Z* after which the integrals tend to increase 
with growing Z. As the studies show, this turn around point is frequently connected 
with appropriate phase functions. This essentially motivates the stopping rule. A 
numerical example is given in Table 1. 

We still mention that for continuous kernels k the bounded linear operators of 
the form F'(x ) : 1^.(0, 1) ->■ L 2 (0, 2) in formula (j2~Tj) denote Frechet derivative of 
the operator (116p at the point xq E L^,(0, 1), which can simply be verified as 

s 

[F'(x )h](s) = J(k(s,q)+k(s,s-q))x (s-q)h(q)dq, 0<s<2, ZiGL 2 (0,l). 



Solving (|16p on a computer requires discretization of the problem, which we did 
as follows. The function values of x are to be reconstructed at N supporting points 
q n which are chosen equidistantly in an interval [<Z m i n , <?max]- The discrete signal is 
denoted by 

x = (xX=i = (*(«0)£U = (A(q n )e^^ =1 (22) 

for n = 1 . . .N and q n = q m \ n + (n — l)Ag with Aq = 9m ^~^ mi " . 
The notation for the output signal is analogous, 

y = (y m )T=-i l = (vism))™-, 1 = (flk^fi 1 (23) 

for m = 1 . . . 2N — 1 and s m = 2g m ; n — q cw + (m — l)Ag. 

Here q cw is the frequency of the quasi-continuous wave ancilla pulse. The kernel 
takes the form 

K=k m)n = K(s m ,q n ,q cw ), m = 1, 2, . . . , 2N - 1, n = 1, 2, . . . , N, (24) 

with K(-,-,-) from (I12p . The autoconvolution operator F is discretized using the 
rectangular rule. That way, 

N 

y m = '^2k(s m ,q j ,q cw )x(q j )x(s m + q cw - qj)Aq, m = 1, 2, . . . , 2N - 1. (25) 

Because of the finite support of x, x(s m + q cw — qj) = for s m + g cw — <fy- < q w n and 
s m + q cw — qj > q mg ,x- The complete operator can be written as a multiplication of 
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a matrix F_(x) 6 C 2N ~ lxN with the vector x, 

y = E.(m)z 

where 



(26) 



F(x) = Aq 



^2,1^2 ^2,2^1 



kN-l,l%N-l fejV-l,2^JV-2 
k N) ix N k N>2 XN-i 
fcjv+l,ia:jv 





V o 













kN,N-l%2 k NtN Xi 
kN+l,N-l%3 kN+l,N-lX2 

k2N-2,N-\XN k2N-2,NXN-l 
k2N-l,NXN ) 



(27) 



Analogously the Frechet derivative F'(xq) is discretized. The m-th entry (F_(xo)h) r 
then reads as 



N 



(F'(x )h) m = J>( )+k( + - qj , Qcw) )x (s m +q cw - qj )h(qj ) Aq. 

3=0 

(28) 

Discretizing the operator L = as 



L 



Aq" 2 



/ 2 -1 ... 0\ 







1 2 -1 ... 
-12 



(29) 



the iteration rule (12111 reads as 



ki+ d := 4) + t (4n)*£ (4)) + £ (4))*(r - ^(4)))- ( 3 °) 



1(0 J 



As initial value x s ^ we take the measured absolute values and a zero phase. 



6. Case studies 



Before turning to examples for the reconstruction of a phase, we want to give an idea 
about the dependency of the measurements on the unknown phase. An example 
is given in Figure 2. There, two slightly different phases are shown together with 
their autoconvolution signal split into phase and absolute values. Both convolutions 
have been performed using the same function for the absolute values A and the 
physical kernel. Although the resulting phases look quite alike, there is a much 
more obvious difference in the resulting absolute values B. The example indicates 
that the complete autoconvolution signal is necessary for the reconstruction. In 
the following we give two examples for the reconstruction. For a more detailed 
discussion about the numerical results and effects we refer to (6pj |. To test our 
reconstruction algorithm, we chose functions for the absolute values A(q) and for 
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Phase I 

4 1 1 1 1 
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frequency (THz) 
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100 200 300 400 500 600 700 800 
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250 300 350 400 450 500 550 600 

frequency (THz) 



convolved phase II 
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frequency (THz) 



convolved absolute values II 



0.8 
0.6 
0.4 
0.2 





100 200 300 400 500 600 700 800 

frequency (THz) 



Figure 2.: In the top row, two different, yet similar looking examples of phases 
for the unknown x are shown. Both are convolved with the same function for the 
absolute values. In the second row, the phases of the convolved pulse y are plotted, 
again being quite alike. However, the absolute values of y, pictured in the bottom 
row, show larger distinctions. 



the phase ip{q). Since the real magnitude of A(q) is unknown, we did most of the 
testing with a maximum value of A max = 10~ 7 so that both x and y were neither 
too large nor too small. Since the kernel adds a mae nitude of 10 28 we arrived 
at approximately 10 14 for the SD-Spider signal. To get reasonable values for the 
regularization parameter, we rescaled it. With a and Aq from the previous section 
we define a = a • • Ag -4 • (10 28 ) -2 . The normalized parameter a will be 

the one given in the evaluation of the experiments. In order to follow the physical 
background of the problem, huge numbers appear inevitable. Both pulses were 
sampled on a fine grid and according to (|26]1 . From this representation we extracted 
absolute values B(s) and phase ip(s), which we again sampled on a coarser grid to 
avoid inverse crime, i.e., an overly good reconstruction that arises as an artifact 
because forward and backward calculations were either performed on the same grid 
or because one grid was a multiple of the other. To all data that is assumed to 
be given as measurements, i.e. the absolute values of both the spectrum of the 
pulse to be reconstructed A(q) and of SD-Spider signal B(s), as well as the phase 
of the SD-signal tp(s), we added to each data point normally distributed random 
noise with zero mean and a standard deviation of 5 percent relative to the correct 
value. The final regularization parameter was chosen according to (|20p . using the 
noisy data A as reference. This means that we assume the solution to be best if 
the deviation of the corresponding absolute values is minimal. 
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Figure 3.: First example using a very smooth fundamental pulse 



6.1. Case I: A very smooth pulse 

As a first example we consider a very smooth pulse shown in Figure 3(a). The 
absolute values possess only one peak and the values of the phase increase from — tt 
and 7r. It should be mentioned that both functions the phase and its first derivative 
are zero in the middle of the frequency domain. Figure 3(b) shows the result of 
the autoconvolution including the added noise. The noise level here is 5 = 5%. The 
best reconstruction was achieved for a = 5.86 • 10 6 . Since we are looking for a very 
smooth solution and use a smoothness penalty for the regularization, the parameter 
is very large. The result is given in Figure 3(c). We were able to recover the phase in 
an acceptable way. Although there are problems for lowest and highest frequencies, 
the middle part is retrieved close to the original phase. Due to the fairly large value 
of the regularization parameter, only few remainders of the highly oscillating noise 
(cf. Figure 3(b)) are still visible. 
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frequent (THi 



(a) Fundamental pulse 





(b) Result of the autoconvolution of the fundamental pulse 





(c) Reconstruction of the fundamental pulse and comparison with the true solution, noise-free, 
a = 2.17 




250 3D0 




(d) Reconstruction of the fundamental pulse and comparison with the true solution with noise 
S = 1%, & = 2.82 



Figure 4.: Second example with a more complicated fundamental pulse. See Table 1 
for characteristic values of certain iterations. 



6.2. Case II: An oscillating phase 

As a second example we have a more complex situation in mind. Here we choose an 
amplitude function with two peaks leading to the noisy function A(q). The phase 
function shown in Figure 4(a) has more of a sinusoidal structure. Figure 4(b) again 
shows the result after the autoconvolution, this time without any noise. Even in this 
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Iteration 


Residual 


Solution smoothness 


Deviation oi absolute values 


(0 


UE.(s.( l) )-y s U 


ll^fi)!! 


II l4,l-Al 


1 


9.5819e-01 


1.6806 


0.5252 


20 


2.4ll5e-02 


4.3935 


0.7916 


40 


2.0682e-02 


4.9909 


0.7937 


60 


l.5369e-02 


4.6496 


0.6077 


80 


2.2482e-03 


3.9273 


0.2563 


100 


1.3792e-03 


3.7857 


0.1964 


120 


1.1022e-03 


3.9429 


0.1701 


140 


9.4595e-04 


4.3199 


0.1623 


142 


9.3083e-04 


4.3642 


0.1622 


143 


9.2340e-04 


4.3865 


0.1622 


144 


o i Rnfio n/i 
y.iDUoe-U4 


4.4uyu 


U.lDZo 


150 


8.7480e-04 


4.5451 


0.1632 


200 


5.7859e-04 


5.8612 


0.1832 


250 


3.1613e-04 


7.2540 


0.2020 



Table 1.: Characteristic values of the iteration process for the solution in Fig. 3(d). 
All values have been normalized to simplify comparison. Since the iteration is 
aimed at minimizing the residuals, the values of the residuals decrease steadily 
during the process. The iteration must be terminated early enough, otherwise the 
solution smoothness gets lost for late iterates. Moreover, the solutions increasingly 
oscillate as a consequence of the ill-posedness. On the other hand, the deviation 
of the absolute values decreases after some starting phase until they reach some 
turning point. At iteration I* = 143 the minimum is reached, and according to our 
stopping rule we choose this as the final solution. If the iteration is continued the 
difference in absolute values increases again, and the reconstructions tend to the 
worse. 

case we need an appropriate parameter a to recover the phase in an acceptable way. 
The noise- free reconstruction is shown in Figure 4(c) with a = 2.17. While in the 
middle part both phases coincide nearly perfectly, problems arise at the boundaries 
as we already observed in the previous example. Because of the structure of the 
autoconvolution equation, the amount of information on the boundaries is much 
lower than in the middle, c.f. (]26p and (]27p . The higher the noise level <5, the 
more severe those problems become. In the last figure a reconstruction for 5 = 1% 
and a = 2.82 is shown where a lot of noise artifacts remained. However, it seems 
that the absolute values are influenced more by the measurement errors than the 
phase. Numerical values for certain iterations are shown in Table 1. Especially 
our criterion to stop the iteration for fixed a is clearly visible there. Namely, at 
= 143, the difference of the absolute values of the iterate to the noisy reference 
values reaches its minimum. The solution at this iteration is the one shown in 
Figure 4(d). 



7. Conclusions 

In this paper, we have studied a new type of kernel-based autoconvolution prob- 
lems, for which a stable approximate solution is required for measurements of 
ultrashort laser pulses with the self-diffraction SPIDER method. The problem is 
formulated as a nonlinear integral equation with complex-valued functions over a 
finite real interval on both sides of the equation. With respect to the mathematical 
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model, the novelty of this inverse problem consists in the occurrence of a physically 
motivated complex kernel function and the observability of the amplitude of the 
incident pulse such that the focus of the paper is on phase reconstruction. After 
reviewing recent developments in ultrashort pulse characterization we have intro- 
duced an appropriate physical model (jlip with kernel (|12p and its mathematical 
description as a nonlinear operator equation. Using this abstract setting, which 
supports the mathematical analysis of the problem, we outlined by Example 13.11 
the local ill-posedness of the inverse problem under consideration. That is, for given 
data the solution cannot be approximated arbitrarily precise, even if the noise on 
the data is sufficiently small. Titchmarsh's convolution theorem formulated for our 
model in Lemma 14.11 enabled us to show that the autoconvolution equation has 
two complex conjugate solutions in case of a trivial kernel. However, the situation 
for arbitrary kernels seems to be an open question. A main goal of this article was 
to suggest and test an adapted regularization approach. We used a variant of the 
Levenberg-Marquardt algorithm pip in the discretized form (|30p . A crucial point 
in the regularization process is the choice of the regularization parameter a. Since 
the amplitude function can be measured as part of the solution, we could motivate 
a specific parameter choice rule (I20p . Under all regularized solutions, which are 
calculated for varying a > 0, we recommend to take the one that approximates 
the observed absolute values in an optimal way. It is known that iterates for reg- 
ularized solutions to ill-posed problems tend to worsen whenever the iteration is 
not stopped early enough. To find an appropriate stopping rule, we monitored the 
deviation between the measured and calculated absolute values, which reaches a 
minimum during the iteration process and then increases again. This minimum is 
used to stop the iteration for fixed a. Table 1 illustrates such behavior for some 
numerical example. Case studies with two different pulses and different noise levels 
indicate opportunities and limitations of our regularization method. 

From a physicist's point of view, the formalism described here opens a perspec- 
tive for using a wider class of nonlinearities for the SPIDER pulse characterization 
method, making it more universally applicable, x nonlinearities such as the self- 
diffraction process are much more widespread than x^ processes that have been 
used nearly exclusively for SPIDER so far. With the solution of the x^ SPIDER 
autoconvolution problem at hand, the integration of a complete optical pulse char- 
acterization set-up in an integrated optical device becomes feasible. As SPIDER 
is the only method that does not require mechanical variation of an optical de- 
lay this enables, in principle, the construction of an optical oscilloscope on chip, 
of course with the requirement to externally process the measured data by the 
discussed regularization procedure. Moreover, there appear other interesting appli- 
cations of x^ SPIDER variants for the characterization of broadband ultraviolet 
pulses, where x^ processes do not offer viable options. For all these intriguing 
applications, we now demonstrated a viable way of retrieving the relevant phase 
data by regularization of the respective autoconvolution problem. 
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